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I.  INTRODUCTION 

An  electrical  discharge  is  a  complex  phenomenon  involving  electrodynamics,  hydrody¬ 
namics,  radiation  processes,  and  chemistry  [Tj.  Over  the  years,  substantial  experimental, 
theoretical,  and  numerical  research  has  been  carried  out  in  this  area  due  to  interest  in  natu¬ 
ral  lightning,  plasma  sources,  and  more  recently,  electrical  discharges  guided  by  femtosecond 
lasers  Em  ME].  Much  of  the  numerical  work  has  focused  on  streamer  propagation,  de¬ 
fined  as  the  propagation  of  a  filamentary  ionization  wave  in  a  cold  gas 

The  approach  typically  used  is  to  couple  Poisson’s  equation  to  a  continuity  equation  for  elec¬ 
tron,  positive  ion,  and  negative  ion  densities,  and  to  solve  these  equations  in  two  dimensional 
axisymmetric  geometry.  Air  chemistry  is  usually  reduced  to  three  reactions:  avalanche  ion¬ 
ization,  attachment,  and  recombination.  A  key  assumption  that  is  almost  universally  made 
is  that  the  rate  coefficients  depend  only  on  the  parameter  E/ng ,  where  E  is  the  electric  held 
and  nfJ  is  the  neutral  gas  density.  Although  this  is  valid  for  a  cold  gas,  it  precludes  accurately 
describing  the  leader  phase  in  which  the  gas  is  hot  enough  to  maintain  the  electron  density 
even  in  low  held  regions.  The  code  described  in  this  report,  SPARC  (Streamer  Propagation 
and  ARCing),  uses  a  more  complete  air  chemistry  model  which  is  appropriate  for  either  the 
streamer  or  leader  phase.  The  code  also  includes  hydrodynamics  of  the  heavy  species  which 
might  be  important  for  long  propagation  distances  in  which  there  is  enough  time  for  the  hot 
gas  to  expand. 

Computer  modeling  of  electrical  discharges  is  difficult  not  only  because  of  the  complex 
physics  involved,  but  also  because  of  the  problem  of  carrying  out  the  calculation  on  a  practi¬ 
cal  timescale.  Streamers  and  leaders  are  characterized  by  a  very  short  spatial  scale  (microns) 
at  their  extremity  and  a  much  longer  spatial  scale  (meters)  in  their  body.  A  similar  disparity 
can  exist  in  the  temporal  scales.  For  example,  in  the  case  of  a  laser  induced  discharge  in 
air,  the  timescale  for  electrons  to  lose  energy  to  Nitrogen  vibrations  is  picoseconds,  while 
the  hydrodynamic  timescale  can  be  microseconds.  For  uniformly  spaced  grid  cells  and  time 
levels,  this  implies  that  millions  of  cells  and  millions  of  time  levels  would  be  needed  just  to 
do  a  one- dimensional  calculation.  The  solution  adopted  in  this  work  is  three  fold.  First,  an 
adaptive  grid  is  used  so  that  tightly  packed  grid  cells  are  only  used  where  they  are  needed. 
Second,  the  code  is  written  to  take  advantage  of  massively  parallel  computer  architectures. 
Third,  a  hcierarchy  of  models  is  introduced  which  allow  various  approximations  to  be  in- 
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voked  optionally  at  run-time.  These  approximations  are  designed  to  allow  longer  time  steps 
to  be  taken  while  maintaining  numerical  stability. 

The  heierarchy  of  models  contains  several  levels  of  approximation.  At  the  highest  level, 
the  full  hydrodynamic  equations  are  solved  for  every  species  including  the  electrons.  This 
requires  resolving  the  plasma  period  and  the  electron  collision  time.  At  the  next  level 
of  approximation,  mobility  limited  flow  is  assumed.  This  requires  resolving  the  charge 
redistribution  time,  but  not  the  plasma  period  or  electron  collision  time.  At  the  next  level 
of  approximation,  the  electron  motion  is  treated  implicitly  and  quasineutrality  is  invoked  in 
certain  evaluations.  In  this  case  none  of  the  detailed  electron  motion  needs  to  be  resolved. 
Instead,  the  minimum  time  step  is  determined  by  fast  chemistry  such  as  electron  cooling. 
The  fast  chemistry  can  also  be  treated  implicitly  if  desired.  In  every  case,  the  heavy  particles 
can  be  held  fixed  for  further  computational  savings. 


II.  FULL  ELECTROSTATIC-HYDRODYNAMIC  MODEL 

A.  Basic  Equations 

For  each  species,  define  the  density  ns ,  fluid  velocity  vs,  energy  density  us ,  mass  ms,  and 
charge  qs.  It  is  assumed  the  translational  and  rotational  modes  are  in  thermal  equilibrium, 
but  the  vibrational  modes  are  not.  The  energy  density  is  therefore  written  as 

us  =  1 nsmsv2s  +  0S  +  (1) 

where  0S  accounts  for  translational  and  rotational  thermal  energy,  and  Ss  accounts  for 
vibrational  energy.  Neglecting  viscosity  and  the  magnetic  field,  the  evolution  of  each  species 


is  described  by  following  coupled  continuity  equations:  m 

dtns  +  V  ■  (nsvs)  =  Ns  (2) 

dt(nsvs)  +  V  ■  (nsvsvs)  =  1E  +  Vs  (3) 

ms  ms 

dtes  +  V  •  (0svs)  =  -Ps V  •  vs  -  V  •  hs  +  Qs  (4) 

dtEs  +  V  •  (S,vs)  =  As  (5) 
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Here,  Ps  is  the  partial  pressure  of  species  s,  E  is  the  electrostatic  field,  and  hs  is  the  heat 
flux.  The  electrostatic  field  is  computed  from 


V20  =  -47 T  J2  Vj 


rii 


E  =  -V0 


The  pressure  is  found  from 


Ps  =  nsTs 

where  the  temperature  Ts  is  related  to  Qs  by 

T  20s 

(3  +  Rs)ns 


(6) 

(7) 

(8) 

(9) 


Here,  Rs  is  the  number  of  rotational  degrees  of  freedom.  The  terms  Ns,  Vs,  Qs,  and  Xs 
result  from  chemical  reactions  and  collisions.  The  part  due  to  reactions  is  denoted  by  the 
same  symbol  primed,  while  the  part  clue  to  collisions  is  denoted  by  a  double-prime.  For 
example, 

Ns  =  N's  +  N”  (10) 


Ohmic  heating,  which  plays  a  crucial  role  in  electrical  discharges,  is  contained  in  Q"s. 


B.  Chemical  Reactions 

Given  a  list  of  chemical  reactions  indexed  by  the  variable  i,  the  reaction  rates  are 

=  (11) 

3 

where  is  a  temperature  dependent  rate  coefficient,  j  varies  over  the  reagents  for  reaction 
i,  and  ri{j)  is  the  stoichiometric  coefficient  for  reagent  j  of  reaction  i.  The  rates  at  are 
assumed  to  be  of  the  form 

OLi  =  ai0 T/8  exp_7i/Ts  (12) 

where  Ts  can  be  the  temperature  of  any  species  involved  in  the  reaction,  and  ai0,  f3i:  and  y* 
are  constants.  The  source  terms  due  to  chemical  reactions  are 


K  =  5Zfti[Pi(a)  _ri(s)' 


(13) 
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v;  =  X!  n<  riU)fisjmjVj  -  ri(s)m8vt 


Qs  =  J2  n>  J2  ri(j)fisj  (—  +  ^s  — )  -  n(s)—  +  €i(s) 


x's=(i-  xs)  -  ^(^)- 


where  pj(s)  is  the  stoichiometric  coefficient  for  product  s  of  reaction  i,  Cj(s )  is  the  partial 
heat  of  reaction  for  species  s,  and  fisj  is  the  fraction  of  the  energy  lost  by  reagent  j  that 
goes  into  species  s.  The  operator  Xs  simply  evaluates  to  zero  for  species  with  vibrational 
degrees  of  freedom,  and  unity  for  species  without  vibrational  degrees  of  freedom.  Note  that 
conservation  of  energy  and  momentum  requires  that  J2S  fisj  —  1- 

As  an  example,  suppose  reaction  i  —  1  corresponds  to  three  body  attachment,  which  has 
the  equation  e~  +  02  +  O2  — 1 >  O 2  +  O2.  Let  the  e~  index  be  j  =  1,  the  02  index  be  j  =  2, 
and  the  O2  index  be  j  =  3.  Then  rqf 1)  =  1  and  ri( 2)  =  2,  so  that  1Z\  =  OL\n\n\.  The 
attachment  energy  is  accounted  for  by  ei(3)  =  0.44  eV.  The  energy  and  momentum  lost  by 
e~  and  02  is  gained  by  (Tj,  so  that  /131  =  /i32  =  1. 

C.  Collisions 

The  part  of  the  source  term  for  the  momentum  equation  due  to  collisions  is 


V"  =  ^ 2njnscrSj  (—  +  — 

\m,  m, 


i  Vj  -  v. 


where  aSj  is  the  cross  section  for  collisions  of  species  s  with  species  j.  For  electron-neutral 
collisions,  the  cross  section  is  taken  to  be  constant  at  5  x  10~15  cm2.  Momentum  is  conserved 
if  the  neutral-electron  cross  section  is  smaller  by  the  mass  ratio.  For  Coulomb  collisions, 


°8j 


4V^SV  In  A 


±  n  u  +  t ; 

ms  rrij  J  \ms  rrij 


where  In  A  is  the  Coulomb  logarithm.  The  heating  term  due  to  collisions  is 


Q  s  f  ^  .i 


pj  -  T.)  +  Y,  ™.v"  •  (v*  V.)  -  •£  x; 


where  vSj  is  the  coefficient  of  thermal  equilibration,  the  second  term  accounts  for  heating 
due  to  friction  (including  ohmic  heating),  and  the  index  /  varies  over  species  that  exchange 
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vibrational  energy  with  the  thermal  energy  of  species  s.  The  thermal  equilibration  coefficient 
is  related  to  the  collision  cross  section  by 


Vsj 


3crsjms 


+ 


Ta 


1/2 


ffl; 


(20) 


evv  evv 


exp 


T 
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T, 


(21) 


rns  +  rrij  \  rns 

The  vibrational  source  term  is 

X"  =  nsev  (l  -  e~£v/q,s)  XuUjV 

V 

where  is  the  vibrational  temperature,  v  is  the  vibrational  level,  ev  is  the  energy  between 
adjacent  vibrational  levels,  j  indexes  the  species  that  excites  the  vibration,  and  Xv  is  a 
two-body  rate  coefficient  which  depends  on  Tj  through 

Xu  =  P  eXP  (“C5 uTj)  (22) 


This  form  was  found  to  give  satisfactory  fits  to  the  tabulated  data  of  Ref.  Ha.  The  oper¬ 
ator  P  multiplies  by  zero  if  the  operand  is  negative  and  by  unity  otherwise  (the  curve  fit 
sometimes  gives  negative  values  for  Tj  0.2  eV).  The  vibrational  temperature  is  found 
from 


vly  — _ v _ 

In  (1  +  evns/Es) 


(23) 


III.  IMPLICIT  ELECTRON  MODEL 


The  full  electrostatic-hydrodynamic  model  discussed  in  the  last  section  requires  that  the 
electron  motion  be  explicitly  resolved.  This  is  computationally  demanding  since  it  requires 
that  the  timestep  be  small  enough  to  resolve  both  the  electron  collision  frequency  and  the 
electron  plasma  frequency.  For  example,  to  accurately  resolve  electron-neutral  collisions 
requires  a  timestep  of  <  100  femtoseconds.  In  this  section  we  give  a  model  that  treats  the 
heavy  species  as  in  the  last  section,  but  uses  more  efficient  equations  for  the  electrons. 

The  need  to  resolve  the  electron- neutral  collision  frequency  can  be  removed  by  replacing 
the  electron  momentum  equation  with  its  equilibrium  solution 

ve  =  qeUeE  ~  VPe  (24) 

vemene 

where  the  subscript  e  refers  to  electrons,  convection  of  momentum  was  neglected,  ve  is  the 
aggregate  collision  frequency  for  electrons  with  all  heavy  particles,  and  we  took 

Ve  =  -veneve 


(25) 
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Unfortunately,  this  does  not  necessarily  improve  the  situation  since  the  charge  redistribution 
time  mei,e/4:7rneql  must  still  be  resolved.  For  electron  densities  higher  than  1015  cm”3  this 
timescale  is  even  shorter  than  the  plasma  period  or  collision  time. 


The  need  to  resolve  the  charge  redistribution  time  emerges  when  Eq.  (24),  the  Poisson 


equation,  and  the  continuity  equation  are  combined.  Let  p  be  the  charge  density  and  j  be 


(26) 


the  current  density.  Inserting  Eq.  (24)  into  j  ~  qeneve  gives 

V  QeUe  J 

where  a  =  qlne/mei'e  is  the  conductivity.  Inserting  this  into  the  equation  of  charge  conser¬ 
vation  dtp  +  V  •  j  =  0  and  using  Poisson’s  equation  V20  =  — 47rp  gives 


+  V0i  •  Ver  =  Eq  •  V <7  —  V  ■  (peVP) 


(27) 


where  0 i  is  the  space  charge  potential.  E0  is  the  known  external  field,  and  pe  =  qe/meve  is 
the  electron  mobility.  The  operator  dt/^+cr  dictates  that  in  regions  of  uniform  conductivity 
relaxation  occurs  on  a  l/47rcr  timescale. 

Conceptually,  the  simplest  way  to  avoid  resolving  the  charge  redistribution  time  is  to 


simply  drop  the  time  derivative  from  Eq.  (27)  and  solve  for  the  equilibrium  0!  every  timestep. 


This  approach  turns  out  to  be  ineffective,  however,  because  the  global  equilibration  rate  is 
not  always  fast  compared  to  other  rates  of  interest.  Another  approach  jTF]  is  to  solve 


Eq.  (27)  using  an  implicit  differencing  scheme  (this  drives  the  equation  toward  its  equilibrium 
solution  no  matter  how  large  the  timestep).  This  results  in  an  elliptical  equation  for  0i 
with  coefficients  that  not  only  vary  in  time,  but  also  have  a  spatial  variation  that  is  not 
generally  separable.  This  is  undesirable  because  the  direct  elliptical  solver  to  be  described 
below  cannot  be  used  on  such  equations.  Furthermore,  it  is  preferable  from  a  programming 
point  of  view  to  use  the  same  field  solver  for  both  the  implicit  electron  model  and  the  full 
electrostatic-hydrodynamic  model. 

Accordingly,  the  approach  adopted  here  is  to  obtain  0  from  the  usual  Poisson  equation, 


but  obtain  p  from  Eq.  (27)  re-written  as  follows: 


(■ dt  +  47ra)  p  =  — E  •  Vct  +  V  •  (/ieVP) 


(28) 


This  equation  can  be  solved  stably  for  timesteps  l/47rcr  by  treating  the  left  hand  side 
implicitly  and  the  right  hand  side  explicitly.  This  results  in  an  equation  that  is  very  easy  to 
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solve.  Further  simplification  can  be  obtained  by  replacing  the  continuity  equation  (|2])  with 
the  quasineutrality  condition 

qene  =  ~  ^ns  (29) 

s^e 

and  the  energy  equation  Q  with 

dtee  =  Qe  (30) 

Note  that  these  replacement  equations  apply  only  to  the  electrons.  Heavy  particle  motion 
is  still  computed  using  Eqs.  (§-<§- 

IV.  NUMERICAL  SOLUTION 

A.  Grid  Convention 

Consider  a  two-dimensional  grid  with  cells  indexed  by  i  and  j,  and  coordinate  axes 
denoted  by  z  and  r.  The  interior  grid  cells  are  labeled  by  the  set  i  G  (1,  2, ...,  Nz)  x  j  e 
(1,  2, ...,  Nr).  The  interior  is  surrounded  by  a  single  layer  of  ghost  cells.  All  quantities 
are  considered  known  at  the  cell  centers.  For  example,  refers  to  the  space  charge 

potential  at  the  center  of  cell  (i,j). 

The  particular  geometry  of  the  grid  does  not  have  to  be  specified  in  formulating  the 
problem.  Instead,  finite  difference  equations  can  be  derived  that  depend  on  the  volumes  of 
the  cells  and  the  areas  of  the  cell  walls.  These  areas  and  volumes  can  be  specified  when  the 
code  is  started,  and  can  even  be  changed  while  the  code  is  running.  This  allows  a  single 
algorithm  to  work  on  either  a  cartesian  or  cylindrical  grid,  and  also  allows  the  grid  to  change 
as  the  simulation  runs.  The  volume  of  cell  (i,j)  is  written  V(i,j).  The  area  of  the  cell  wall 
bounding  cell  (i,j)  on  the  negative  z-side  is  Az(i,j).  The  area  of  the  cell  wall  bounding  cell 
(i,j)  on  the  negative  r-side  is  Ar(i,j).  The  position  of  the  center  of  cell  (i,j)  is  (zi,rj).  It  is 
assumed  that  the  i-position  of  any  cell  is  independent  of  j,  and  vice-versa.  This  is  equivalent 
to  the  statement  that  the  grid  lines  must  be  orthogonal. 


B.  Second  Order  Time  Advance 


The  differential  equations  describing  the  evolution  of  all  the  quantities  considered  in  the 
SPARC  model  can  be  put  in  the  form 


af  =  Tf 


(31) 


where  f  is  the  vector  of  all  quantities  to  be  advanced  and  T  is  an  operator  that  does  not 
depend  explicitly  on  time.  Let  the  value  of  f  at  time  level  n  be  denoted  fn.  To  compute 
fn+1  with  second  order  accuracy,  a  variant  of  the  midpoint  method  is  used.  In  the  usual 
midpoint  method,  a  provisional  value  fn+1/2  is  computed  using 


|"n+ 1/2  jra  I  ^  7 

~2 

and  the  value  at  the  new  time  level  is  computed  using 

pi+l  _  _|_  AtTfn~*~1/2 


(32) 


(33) 


This  description  is  not  quite  adequate  for  SPARC  because  of  the  subtleties  involved  in 
evaluating  certain  equations.  However,  the  midpoint  method  can  be  easily  generalized  as 
follows.  Suppose  we  have  an  operator  F(At,  g)  that  advances  fn  by  a  time  interval  At  using 


g  as  the  vector  of  quantities  to  be  used  in  evaluating  the  right  hand  side  of  Eq.  (31 ).  Suppose 
further  that  if  g  =  fn,  the  advance  is  first  order  accurate: 

fn+1  =  F(At,  fn)f"  +  O(At)  (34) 

Then,  second  order  accuracy  is  achieved  by  using 

P+1/2  =  fn)  in  +  O(At)  (35) 

fn+i  =  fn+ 1/2)P  +  o(At2)  (36) 


It  remains  only  to  identify  the  components  of  F  with  the  first  order  algorithms  that  are  used 
to  advance  the  quantities  of  interest. 


C.  Finite  Difference  Form  of  Divergence  and  Laplacian 

The  finite  difference  form  of  the  divergence  of  a  discrete  vector  field  E (i,j)  defined  on 
any  grid  where  the  geometry  is  specified  as  above  is 

V  ■  E  =  —  DiEz(i  —  1/2,  j)  +  +  1/2,  j)  —  D3Er(i,j  —  1/2)  +  D^Er[i:j  +  1/2)  (37) 
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where  the  subscripts  on  E  select  a  vector  component,  and 


=  Az(i,j)/V(i,j)  (38) 

D2(i,j)  =  Az(i  +  1  ,j)/V(i,j)  (39) 

Ds(i,j)  =  Ar(i,j)/V(i,j)  (40) 

D,i(i,j)  =  +  1  )/V{i,j)  (41) 


Note  that  the  arguments  of  Di . . .  Z/4  were  suppressed  in  Eq.  (37).  Also,  fractional  indices 
indicate  that  an  average  is  to  be  taken  over  the  values  corresponding  to  the  two  nearest 
integral  indices.  For  example, 


Ez(i  -  1/2,  j)  =  E,(i  -  1, j)/2  +  Ez(i,j)/ 2 


(42) 


The  form  of  the  Laplacian  consistent  with  this  definition  of  the  divergence  is  found  by 
substituting 

Ez(i  -  1/2,  j)  =  (0(f  -  1  ,j)  -  (j)(i,j))/Ai_  1/2  (43) 

Ez{i  +  1/2,  j)  =  ((j>(i,j)  -  4>(i  +  1,  j))/Ai+1/2  (44) 

Er{i,j  -  1/2)  =  (0(i,  j  -  1)  -  <j>(i,j))/ Aj-i/2  (45) 

Er{i,j  +  1/2)  =  -  <p(i,j  +  l))/Ai+i/2  (46) 

into  the  formula  for  the  divergence.  Here  we  have  defined 


1/2  ^z+1 

Aj+ 1/2  =  rj+ 1  -  rj 


This  gives  for  the  Laplacian 


Dl(j)(i-l,j)  D2(j){i  +  l,j)  D3cj)(i,j- 1)  D4cj)(i,j  + 1) 

V  0  =  - 7 - t - 7 - 1 - 7 - 1 - 7 - b  D5(j)(l,j) 


A,- 


i- 1/2 


A,-. 


*+1/2 


A 


1-1/2 


Aj+1/2 


where 


£5 


-Di  ^  ^2  -D3  \ 

Aj_i/2  Aj+1/2  Aj-l/2  Aj+1/2  / 


(47) 

(48) 


(49) 


(50) 
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D.  Poisson  Solver 


The  Poisson  solver  used  in  SPARC  takes  advantage  of  the  fact  that  for  the  problem  of 
interest,  Nz  Nr,  where  we  now  associate  z  with  the  axial  direction  and  r  with  the  radial. 
Even  with  an  adaptive  grid,  the  number  of  axial  cells  is  still  much  greater  than  the  number  of 
radial  cells  because  the  cell  size  has  to  be  changed  gradually  to  maintain  accuracy.  SPARC 
takes  advantage  of  Nz  3>  Nr  by  solving  the  Poisson  equation  using  a  technique  that  requires 
0(NzN?)  operations.  Since  Nr  is  not  too  large  (~  100),  the  technique  works  well.  In  fact, 
it  turns  out  that  the  Poisson  solver  does  not  dominate  the  computation  time  in  most  cases. 
Unlike  iterative  methods,  the  approach  developed  here  gives  the  solution  to  within  machine 
precision  in  a  fixed  number  of  operations.  It  also  parallelizes  well  via  domain  decomposition 
in  the  axial  direction. 

The  finite  difference  form  of  the  Poisson  equation  can  be  put  in  the  form  of  a  matrix 
equation  if  (f>(i,j)  and  p(i,j )  are  regarded  as  matrices  where  i  is  the  axial  (column)  index 
and  j  is  the  radial  (row)  index.  The  matrix  form  of  the  Poisson  equation  is 

R0+ (Z0T)7  = -4vrp  (51) 


where  R  is  the  radial  part  of  the  Laplacian  operator  written  in  matrix  form  and  Z  is  the 
axial  part.  The  superscript  T  indicates  that  the  transpose  is  to  be  taken.  The  matrices  R 
and  Z  are 

L>i(f)/Ai_i/2  m  =  i-  1 

-Di(i)/ —  D2(i)/Ai+i/2  m  —  i 
D2{i)/Ai+ 1/2  m  =  i  +  1 

0  otherwise 


^im  \ 


(52) 


Rj  rri  \ 


(53) 


D3(j)/Aj_1/2  m  =  j-  1 

-D3(j)/Aj_1/2  -  D4(j)/ Aj+1/2  m  =  j 

D4(j)/Aj+1/2  m  =  j  + 1 

0  otherwise 

Here  the  requirement  emerges  that  D i  and  D2  be  independent  of  j,  and  D3  and  D4  be 
independent  of  i.  This  requirement  will  always  be  met  provided  the  grid  lines  are  orthogonal. 
We  now  arrange  the  eigenvectors  of  R  in  columns  to  form  a  matrix  denoted  by  H This 
matrix  has  the  property  that 

HRH~l  =  A  (54) 
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where  A  is  diagonal.  Multiplying  the  Poisson  equation  by  H  on  the  left  hand  side  causes 
the  problem  to  decouple  into  Nr  independent  one  dimensional  problems: 


D  i 


A,_ 


i- 1/2 


-0(d  -  1  ,j) 


D  i 


A,;_ 


+ 


Do 


i- 1/2 


A 


i+l/2 


At  4>{hi)  + 


Do 


A. 


i+l/2 


-0(*  +  l,j)  =  -4np(i,j)  (55) 


where  <fi  =  H<p  and  p  =  Hp.  This  equation  can  be  solved  by  the  usual  tridiagonal  algorithm. 
This  leads  to  the  following  procedure  for  obtaining  <p  given  p: 


1.  Diagonalize  R  and  store  the  eigenvalues  A j  and  the  matrices  H  and  H~l.  We  do  this 
using  the  QR  algorithm  and  inverse  iteration.  Since  the  radial  grid  spacings  do  not 
evolve,  R  stays  constant  throughout  the  simulation  and  this  can  be  done  once  at  the 
beginning. 


2. 

3. 

4. 


Form  the  transformed  charge  density  p  =  Hp.  This  takes  0(NzN. (?)  operations. 


Solve  the  tridiagonal  system  (55)  for  <p.  This  takes  0(NzNr )  operations. 


Form  the  solution  in  real  space,  <p  —  H  l<p.  This  takes  0(NzNp)  operations 


Note  that  if  the  bounded  discrete  space  considered  here  is  replaced  by  an  unbounded  con¬ 
tinuous  space,  H  becomes  the  Hankel  transform  of  zero  order. 


E.  Fluid  Advance  and  Chemistry 

The  fluid  equations  (|2])-(|5j)  are  advanced  using  the  well  known  method  of  flux  corrected 
transport  (FCT).  The  particular  variant  of  FCT  used  is  essentially  a  C++  rewrite  of  the 
LCPFCT  (NRL  Laboratory  for  Computational  Physics  and  Fluid  Dynamics  FCT)  routine 
im  The  routine  is  extended  to  two  dimensions  using  the  operator  splitting  approach 
described  in  Ref.  mt 

The  source  terms  for  the  fluid  equations  include  terms  following  directly  from  the  Boltz¬ 
mann  equation,  such  as  the  heat  flux  term,  as  well  as  terms  that  are  put  in  “by  hand”  to 
account  for  chemical  reactions.  These  terms  can  all  be  found  by  straightforward  differencing. 
Of  course,  the  differencing  has  to  account  for  variations  in  the  grid  geometry  in  the  manner 
described  above. 
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F.  Sliding  Rezone 


A  streamer  generally  consists  of  an  elongated  region  called  the  “body”  in  which  the  axial 
derivatives  dz  are  small,  and  a  much  shorter  region  at  the  end  of  the  body  (the  “tip” )  where 
dz  is  large.  The  sliding  rezone  technique  minimizes  the  number  of  grid  cells  needed  to  model 
the  streamer  by  using  large  cells  in  the  body  and  small  cells  in  the  tip.  The  technique  works 
by  keeping  the  total  number  of  axial  grid  cells,  Nz,  constant,  while  varying  the  spacing 
between  the  grid  cells  as  the  simulation  runs. 

The  function  giving  the  length  of  the  ith  grid  cell  is 


1  +  Af  (n- l) 

0  <  i  <  N 

“Regionl” 

i  +  MAd) 

N  <i  <2N 

“Region2” 

1 

2N  <  i  <  3N 

“Region3” 

3N  <  i  <  AN 

“Region4” 

where  Az  is  the  minimum  cell  length,  N  =  Nz/ 4  is  the  number  of  cells  in  each  region, 
and  the  function  /  is  a  smooth  polynomial  that  ramps  from  zero  to  one  and  back  to  zero 
according  to 


f(j)  =  { 


10(2r)3  —  15(2r)4  +  6(2r)5 
1  —  10(2r  —  l)3  +  15(2t  —  l)4 
0 


6(2r 


!)5 


0  <  t  <  1/2 
1/2  <  r  <  1 
otherwise 


The  grid  is  modified  as  the  simulation  runs  by  varying  the  parameters  A,  B ,  and  C.  Thus, 
the  resolution  in  regions  1,  2,  and  4  can  be  changed,  but  region  3  is  always  kept  at  the 
highest  resolution. 

Region  1  is  used  to  contain  the  initial  streamer.  It  does  not  evolve  during  the  course 
of  the  simulation.  It  contains  a  high  resolution  region  at  both  ends.  The  high  resolution 
sub-region  near  the  electrode  (i  —  1)  is  not  currently  needed,  but  can  be  used  in  the  future 
to  support  a  sheath  model.  The  sub-region  near  i  =  N  is  kept  at  high  resolution  due  to  the 
fact  the  streamer  tip  leaves  perturbations  at  its  initial  location  long  after  it  has  propagated 
away.  The  constant  characterizing  region  1,  A,  is  specified  at  the  beginning  of  the  simulation 
and  does  not  have  to  be  changed  thereafter.  In  particular, 


A 


1 

S 


(56) 
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where 


N 


s  =  Y,f 


i=  1 


i  -  1 
AT-  1 


(57) 


and  L\  is  the  desired  length  of  region  1. 

Region  2  is  used  to  contain  the  streamer  body.  Early  in  the  simulation,  the  tip  propagates 
through  this  region,  so  the  highest  resolution  is  needed.  Therefore  the  parameter  B  is  initially 
set  to  zero.  Once  the  tip  reaches  the  center  of  region  3,  however,  B  is  increased  such  that 
the  tip  stays  centered  in  region  3.  In  particular,  when  the  code  detects  that  the  peak  of  the 
axial  electric  held  has  moved  more  than  one  grid  cell  beyond  the  center  of  region  3,  a  new 
value  of  B  is  calculated  using 


B 


I  (  z° 

S  VAx: 


5  N 

~Y 


+  1  -  as) 


(58) 


where 

x0  =  Az  +  AS  +  B'S^j  (59) 

and  B’  is  the  old  value  of  B.  The  new  value  of  B  is  chosen  to  shift  the  cells  in  region  3  by 
an  amount  Ax:.  The  parameter  Zq  is  the  location  of  the  center  of  region  3  before  the  shift. 

Region  3  stays  at  the  highest  resolution  throughout  the  simulation,  and  is  used  to  contain 
the  streamer  tip.  This  is  the  high  resolution  region  that  “slides”  forward  with  the  tip  as  the 
parameter  B  is  increased. 

Region  4  is  used  to  model  the  space  between  the  streamer  tip  and  the  “ground”  electrode. 
At  the  start  of  the  simulation  this  region  is  at  its  largest,  and  the  parameter  C  is  chosen 
to  make  the  end  of  the  simulation  box  coincide  with  the  desired  position  of  the  ground.  In 
particular, 

C  =l  -4  N-AS-  BS^j  (60) 

where  Ltot  is  the  total  length  of  the  simulation  box.  This  has  to  be  re-evaluated  every  time 
B  changes.  That  is,  in  order  to  keep  Ltot  constant,  C  has  to  be  reduced  every  time  B  is 
increased 


G.  Boundary  Conditions 

For  the  Poisson  equation  and  the  hydrodynamic  equations  boundary  conditions  are 
needed.  For  the  hydrodynamic  quantities  the  normal  derivative  is  taken  to  vanish  at  all 
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boundaries.  For  the  virtual  boundary  at  r  =  0,  the  normal  derivative  vanishes  for  scalars 
and  the  quantity  itself  vanishes  for  vectors.  These  boundary  conditions  allow  fluid  to  flow 
freely  into  or  out  out  of  any  boundary.  In  reality,  there  would  be  a  sheath  near  the  electrode 
which  we  neglect. 

In  the  case  of  the  Poisson  solver,  the  boundary  conditions  are  chosen  so  that  in  the 
absence  of  space  charge,  the  field  is  consistent  with  that  produced  by  a  spherical  electrode. 
However,  this  can  only  be  done  approximately  because  of  the  shape  of  the  grid  lines.  In 
particular,  the  potential  in  any  cell  is 


4>{hj)  =  Mh  j)  +  Mhj)  (61) 

•  .n.  _  Re4>00 

\Jr2  +  ( Re  +  Zi )2 

where  0i  is  the  space  charge  potential,  0 o  is  the  potential  in  the  absence  of  space  charge, 
Re  is  the  radius  of  the  electrode,  and  0O o  is  the  potential  at  the  surface  of  the  electrode. 
The  expression  for  0O  is  taken  to  be  valid  both  in  the  interior  and  on  the  boundaries,  while 
0i  is  assumed  to  vanish  at  the  boundaries.  For  the  virtual  boundary  at  r  =  0,  the  normal 
derivative  vanishes.  This  formulation  is  valid  if  rj  <C  Re  for  all  j  since  then  the  surface  of 
the  electrode  nearly  coincides  with  the  planar  boundary  of  the  grid. 

A  real  air  discharge  can  often  be  considered  radially  unbounded.  It  is  therefore  desirable 
to  move  the  radial  simulation  boundary  as  far  from  the  axis  as  possible.  In  order  to  do 
this  while  keeping  the  number  of  radial  grid  cells  small,  the  cell  spacings  are  increased  with 
increasing  r.  In  particular,  the  radial  cell  size  A  rj  is  taken  to  vary  as 


Ar-j  = 


(63) 


Ar0  0  <j<Nu 

Ar0(l  +  8r)j-Nu  Nu  <  j  <  Nr  +  1 
where  Ar0  determines  the  highest  resolution,  Nu  is  the  number  of  cells  in  the  region  of 
uniform  resolution,  and  8r  is  a  number  which  determines  how  fast  the  cells  grow  for  j  >  Nu. 
Note  that  the  ghost  cells  j  =  0  and  j  =  Nr  +  1  are  included  in  the  definition.  For  the 
examples  in  this  report  we  use  Sr  =  0.1. 


H.  Parallelization 

Parallelization  of  SPARC  is  accomplished  via  two-dimensional  domain  decomposition. 
That  is,  the  two-dimensional  grid  is  divided  into  Np  =  Npz  x  Npr  subgrids.  The  total  number 
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of  subgrids  Np  is  usually  set  equal  the  number  of  available  processors  so  that  each  processor 
calculates  all  the  quantities  on  a  single  subgrid.  On  processors  that  support  simultaneous 
multi-threading  (SMT,  sometimes  called  “hyperthreading”)  it  is  advantageous  to  to  assign 
multiple  subgrids  to  each  processor. 

Since  the  equations  on  each  subgrid  are  not  independent  of  one  another,  information 
has  to  be  exchanged  between  processors.  The  turboWAVE  framework,  on  which  SPARC  is 
built,  supports  two  communications  methods.  For  shared  memory  systems,  POSIX  threads 
(“ptrheads”)  are  used.  For  distributed  memory  systems,  the  Message  Passing  Interface 
(MPI)  is  used. 

The  SPARC  module  most  difficult  to  parallelize  via  domain  decomposition  is  the  Poisson 
solver.  This  is  because  the  potential  in  any  grid  cell  depends  on  the  charge  density  in  every 
other  grid  cell.  In  order  to  divide  the  domain  in  the  radial  direction,  the  well  known  transpose 
technique  is  used.  To  divide  the  domain  in  the  axial  direction,  we  use  the  turboWAVE 
parallel  tridiagonal  solver  described  in  Ref.  [IB]. 

The  hydrodynamics  module  parallelizes  via  domain  decomposition  in  a  straightforward 
way.  This  is  because  each  grid  cell  is  updated  using  only  information  from  itself  and  its 
nearest  neighbors  (in  one  evaluation  the  next-nearest  neighbor  is  also  needed).  Hence,  the 
cells  in  the  interior  of  a  subgrid  can  be  advanced  one  time  level  independently  of  the  other 
subgrids.  The  ghost  cells  are  then  obtained  by  receiving  the  corresponding  interior  cell  from 
the  adjacent  subgrid.  This  exchange  of  information  is  carried  out  in  a  separate  sweep  for 
each  direction  of  the  domain  decomposition. 

Finally,  the  chemistry  module  parallelizes  trivially.  That  is,  the  reaction  rate  in  a  cell 
depends  only  on  quantities  in  the  same  cell.  Hence,  each  subgrid  can  update  itself  indepen¬ 
dently  of  the  others. 


V.  BENCHMARKING 

Due  to  the  complexity  of  SPARC,  it  is  important  to  benchmark  as  many  of  the  individual 
computational  elements  of  the  code  as  possible.  In  the  following  sections  we  show  that  the 
individual  components  give  the  expected  results  in  cases  where  either  an  analytical  result  is 
known,  a  general  trend  is  expected,  or  a  comparison  with  another  code  can  be  made. 
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FIG.  1:  Comparison  of  simulated  and  theoretical  radial  electric  fields  due  to  a  long  cylinder. 

A.  Poisson  Solver 

The  simplest  test  to  perform  is  to  verify  that  the  Poisson  solver  gives  the  correct  elec¬ 
trostatic  field  in  some  analytically  tractable  limiting  cases.  As  an  example,  Fig.  [l]  shows 
the  simulated  and  theoretical  radial  electric  field  due  to  a  uniformly  charged  cylinder  with  a 
radius  small  compared  to  its  length.  The  agreement  is  satisfactory.  We  have  also  compared 
the  fields  calculated  by  SPARC  with  those  calculated  using  the  Fourier  Analysis  and  Cyclic 
Reduction  (FACR)  technique  and  found  the  two  results  to  be  nearly  indistinguishable.  We 
further  verified  that  the  results  obtained  on  a  non-uniform  grid  are  consistent  with  those 
obtained  on  a  uniform  grid. 


B.  Chemistry 

Chemical  reactions  in  SPARC  are  defined  in  the  input  file  and  not  by  the  code  itself. 
Specifically,  when  the  code  starts  it  reads  the  input  file  and  creates  an  instance  of  the  same 
C++  object  for  each  reaction  defined  therein.  Chemical  species  are  handled  the  same  way. 
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Hence,  to  benchmark  the  objects  pertaining  to  chemistry  it  is  sufficient  to  check  just  a  few 
reactions.  For  example,  we  checked  that  the  ionization  rate  was  computed  correctly  for 
various  values  of  the  held  and  gas  density.  A  more  interesting  test  was  to  verify  that  the 
electron  cooling  rate  was  consistent  with  that  predicted  by  CHMAIR  II.  In  the  initial  testing 
it  was  assumed  that  nitrogen  vibrations  were  the  main  sink  of  electron  energy.  This  turned 
out  to  be  a  poor  assumption  when  the  electron  temperature  is  high,  as  illustrated  in  Fig.J^a). 
By  accounting  for  the  energy  lost  to  electronic  excitations  into  triplet  states,  the  much 
improved  agreement  shown  in  Fig.  |2](b)  was  obtained.  Most  of  the  remaining  disagreement 
appears  to  be  due  to  the  use  of  different  coefficients  for  the  vibrational  excitations.  SPARC 
currently  uses  the  tables  from  Ref.  |T5] ,  while  CHMAIR  II  uses  those  from  Ref.  [19]. 

C.  Electrodynamics 

The  primary  electrodynamic  effect  relevant  to  streamer  and  leader  propagation  is  the  de¬ 
velopment  of  a  charged  tip  which  leads  to  held  enhancement.  To  observe  this  held  enhance¬ 
ment  in  SPARC,  we  initialize  a  plasma  filament  into  the  simulation  box  with  all  chemical 
reactions  turned  oh,  and  apply  a  constant  external  electric  held.  The  density  in  the  plasma 
filament  is  of  the  form  ne(z,r)  =  nsf  (z)g{r)  where 


1 

0  <  z  <  Ls 

m  = 

1  - 

(IOC3  -  15C4  +  6C5) 

Ls  <  z  <  Lt 

0 

/-  z  ~  Ls 

Lt~Ls 

otherwise 

1 

0  <  r  <  Rs 

g(r)  =  < 

1  - 

(IOC3  -  15C4  +  6C5) 

Rs  <  r  <  Rt 

0 

CO 

a? 

I  ' 

oC 

otherwise 

This  represents  a  uniform  cylinder  of  length  Ls  and  radius  Rs ,  together  with  axial  and  radial 
transition  regions  of  dimension  Lt  —  Ls  and  Rt~Rs,  respectively.  We  vary  Ls  and  ns,  while 
holding  hxed  Lt  —  Ls  =  25  pm,  Rt  —  Rs  =  25  pm,  and  Rs  =  25  pm.  The  filament  is  initialized 
in  contact  with  the  electrode  where  it  is  assumed  all  normal  derivatives  vanish  (we  ignore 
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FIG.  2:  Comparison  of  electron  cooling  rate  as  computed  by  SPARC  and  CHMAIR  II:  (a)  before 
including  electronic  excitations  (b)  with  electronic  excitations. 
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FIG.  3:  Field  enhancement  factor  for  various  aspect  ratios  and  plasma  densities. 

the  sheath).  The  aspect  ratio  is  defined  as  a  =  Ls/Rs.  Figure [3] plots  the  held  enhancement 
factor  as  a  function  of  the  aspect  ratio  for  three  plasma  densities.  As  expected,  the  held 
enhancement  factor  is  of  the  order  of  the  aspect  ratio.  It  increases  with  both  aspect  ratio 
and  plasma  density. 

Another  electrodynamics  benchmark  that  was  performed  was  to  measure  the  held  re¬ 
laxation  time  in  the  body  of  the  plasma  filament.  We  verified  that  the  held  varies  as 
E  =  where  E0  is  the  initial  applied  held,  and  t0  ~  RC  is  a  constant  that  depends 

on  the  plasma  density  and  geometry.  Here,  RC  is  the  product  of  the  resistance  and  capaci¬ 
tance  of  the  plasma  filament.  The  RC  constant  for  a  large  aspect  ratio  uniformly  conducting 
cylinder  is 

a 2 

RC  «  2 tr- —  (64) 

In  a 

where  tr  =  1/Ana  is  the  charge  redistribution  time  and  a  is  the  aspect  ratio.  Figure  [4]  shows 
that  the  simulated  held  decays  at  this  rate  to  within  a  factor  of  two.  The  factor  of  two 
discrepancy  is  not  surprising  given  the  approximate  nature  of  the  theoretical  expression. 
Figure  [5]  shows  a  similar  comparison  holding  the  aspect  ratio  hxed  at  10  and  varying  the 
electron  density  (and  hence  the  conductivity). 
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FIG.  4:  Field  relaxation  time  vs.  aspect  ratio  for  an  electron  density  of  2.5  x  1016  cm-3.  The 
corresponding  charge  redistribution  time  is  0.012  ps. 


1014  1015  1016  1017 


Electron  Density  (cm  3) 

FIG.  5:  Field  relaxation  time  vs.  electron  density  for  an  aspect  ratio  of  10. 
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D.  Adaptive  Grid 


The  sliding  rezone  algorithm  was  tested  by  performing  a  simulation  that  can  be  done  both 
on  a  uniform  grid  and  on  an  adaptive  grid  and  comparing  the  two  results.  This  simulation 
models  streamer  propagation  under  the  assumption  that  the  ionization  rate  depends  on 
E/rig  where  E  is  the  electric  field  and  ng  is  the  gas  density  (this  assumption  is  removed 
later).  Also,  a  generic  molecule  M2  is  used  to  represent  either  oxygen  or  nitrogen.  The 
parameters  of  the  simulation  are  given  in  Table  [TJ  For  this  simulation  and  for  the  rest  of  this 
report,  the  functional  form  of  the  initial  plasma  density  is 


n, 


:(z,  r)  =  n,f(z)e  r‘,ri  +  nhe  r'/rt  +  nb 


(65) 


where 


1 

0  <  z  <  Ls 

f(z)  =  < 

^  rts—nf 

ns 

(IOC3  -  15C4  +  6C5) 

Ls  <  z  <  L 

k  nf/ns 

z-Ls 

Lt  -  Ls 

otherwise 

These  parameters  are  interpreted  as  follows:  ns  is  the  plasma  density  in  the  streamer,  rif 
is  the  density  of  the  pre-ionization  ahead  of  the  streamer,  is  a  uniform  background,  Ls 
is  the  length  of  the  streamer  body,  Lt  is  the  length  of  the  streamer  body  plus  tip,  and  r0 
is  the  streamer  radius.  The  second  term  can  be  used  to  include  a  halo  due  to  the  radiation 
surrounding  the  streamer,  although  in  this  example  rih  =  0.  The  result  of  the  benchmark  is 
shown  in  Fig.  [6]  The  results  from  the  adaptive  grid  are  nearly  identical  to  those  from  the 
uniform  grid.  However,  the  uniform  grid  used  106  grid  cells  to  model  a  5  cm  long  region, 
while  the  adaptive  grid  used  only  105  grid  cells  to  model  a  6.5  cm  long  region. 


E.  Implicit  Model  vs.  Explicit  Model 

The  final  benchmark  considered  in  this  report  is  a  comparison  of  the  implicit  electron 
model  with  the  full  electrostatic-hydrodynamic  model.  The  two  models  are  compared  in 
terms  of  their  predictions  regarding  streamer  propagation  in  pure  nitrogen.  The  parameters 
of  the  two  simulations  are  shown  in  Table  |TT]  The  simulations  include  heavy  particle  motion. 
Both  calculations  were  run  on  32  processors  of  the  IBM  cluster  1600  “Kraken”  at  the  Naval 
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TABLE  I:  Parameters  for  Adaptive  Grid  Benchmark 


Parameter 
Electric  Field 
Electron  Mobility 
ns,  Streamer  Density 
ro,  Streamer  Radius 
Ls ,  Streamer  Length 
rif,  Preionization  Density 
rih,  Halo  Density 
rib j  Background  Density 
Ionization  Ratea 
Attachment  Rate 
Recombination  Rate 
Nz ,  Axial  Cells 
Nr,  Radial  Cells 
Nu.  Uniform  Cells 


Value 
24  kV / cm 
1800  cm2/V-s 

2.5  x  1016  cm”3 
100  n m 
1  cm 

2.5  x  1015  cm"3 
0 

2.5  x  109  cm-3 
5  x  io-8e_24°/E  cm3/s 
2.2  x  10_3°  cm6/s 
10- 7  cm3/s 

1000 

100 

50 


Comment 


Field  is  uniform 


see  Eq.  ( 65 ) 


1/e  Definition 
Lt  —  Ls  =  100  fj, m 


see  Eq.  ( 65 ) 


see  Eq.  ( 65 ) 


e“  +  M2  ->•  2e“  +  M2+ 
e”  +  2M2  -»•  M-2  +  M2“ 
e~  +  M/~  —■ >  M2 
104  for  uniform  grid 


see  Eq.  ( 63 ) 


aE  is  the  electric  field  in  kV/cm 

Oceanographic  Office.  The  explicit  calculation  ran  for  20  hours  while  the  implicit  calculation 
ran  for  37  minutes.  A  comparison  of  the  electron  temperature  after  21  ns  is  shown  in 
Fig.0  The  streamer  tip  can  be  identified  as  the  hottest  region  on  the  axis.  Evidently, 
the  implicit  algorithm  predicts  streamer  velocities  that  are  in  agreement  with  the  explicit 
algorithm  to  within  about  5%.  An  interesting  feature  of  the  explicit  result  is  the  discontinuity 
that  appears  near  the  radial  cell  index  70  for  z  >  5  cm.  A  study  of  the  time  evolution 
of  this  discontinuity  reveals  that  it  is  a  thermal  shockwave  propagating  radially  outward. 
This  shockwave  can  actually  be  seen  in  the  implicit  calculation  also,  but  it  is  weaker  and 
appears  only  at  earlier  times.  This  suggests  that  the  implicit  algorithm  introduces  numerical 
dissipation  into  the  system. 

The  subtle  differences  between  the  implicitly  and  explicitly  computed  electron  tempera¬ 
ture  could  be  due  to  several  factors.  First,  explicit  differencing  is  generally  more  accurate 
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FIG.  6:  Comparison  of  simulations  using  adaptive  and  uniform  grids.  The  on-axis  electric  field  at 
t  =  800  ps  is  plotted  vs.  the  axial  coordinate,  z,  for  both  cases. 

than  implicit  differencing.  However,  it  was  found  that  repeating  the  implicit  calculation 
with  the  time-step  reduced  by  a  factor  of  four  gave  the  same  result.  Second,  the  explicit 
calculation  includes  convection  of  electron  momentum  and  energy  density.  This  can  occur 
in  the  tip  where  gradients  are  large,  although  simple  estimates  suggest  the  effect  is  not  sig¬ 
nificant.  Third,  the  explicit  calculation  allows  for  plasma  waves  since  it  does  not  assume  the 
momentum  equation  is  in  equilibrium.  Finally,  the  implicit  calculation  uses  the  quasineutral 
approximation  to  evaluate  the  electron  density  for  use  in  computing  chemical  reaction  rates. 
This  could  result  in  large  errors  in  the  reaction  rates  in  a  small  region  near  the  tip. 


VI.  SAMPLE  RESULTS 

In  this  section  we  give  results  from  three  SPARC  simulations  which  use  the  implicit 
electron  model.  In  the  first  example  a  streamer  is  propagated  over  a  distance  of  40  cm  in 
the  absence  of  attachment.  In  the  second,  a  streamer  to  leader  transition  is  observed  with 
attachment  included.  In  the  third,  the  heavy  species  are  allowed  to  move  and  hydrodynamic 
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TABLE  II:  Parameters  for  Implicit  vs.  Explicit  Benchmark 


Parameter 

Electrode  Voltage 

Electrode  Diameter 

Electron-Neutral  Cross  Section 

Ion-Neutral  Cross  Section 

Electron-Ion  Cross  Section 

ns,  Streamer  Density 

ro,  Streamer  Radius 

Ls,  Streamer  Length 

rif,  Preionization  Density 

rih,  Halo  Density 

rh,  Halo  Radius 

rib ,  Background  Density 

Ionization  Ratea 

Excitation  Rate6 

De-excitation  Rate 

Recombination  Rate 

Vibrational  Excitation 

Nz.  Axial  Cells 

Nr,  Radial  Cells 

Nu ,  Uniform  Cells 

At,  Time  Step 


Value 


500  kV 
40  cm 

5  x  10" 15  cm2 
cnP 


2  x  10”17 


see  Eq.  (18) 

1016  cm-3 
40  /im 
1  cm 

10 15  cm-3 
109  cnr3 
1  cm 

200  cm~3 

1.6  x  10~8Te1/,2e_1'’2/Te  cm3/s 

5.4  x  io_7Te_0-32e_9-52/Te  cm3/s 

3  x  10~9  cm3/s 

4.3  x  10~8Te_o'39  cm3/s 

see  Ref.  [15] 

800 

150 

50 

1.42  ps  (implicit),  0.07  ps  (explicit) 


Comment 
constant  voltage 
spherical  electrode 


see  Eq.  (65) 


1/e  Definition 
Lt  —  Ls  =  100  /mi 


see  Eq.  (65) 


see  Eq.  (65) 


1/e  Definition 


e^  +  N2  ->  2e-  +  JV2+ 
e“  +  N2  — >  e~  +  JV2[*] 
[*]  +  N2  — >  2N2 
e~  +  V+  N2 


see  Eq.  (63) 


“Temperatures  are  in  eV 
hV2[*]  represents  any  triplet  state 

expansion  is  observed.  I11  all  the  examples  in  this  report,  the  initial  temperature  is  0.025  eV 
for  heavy  species,  and  1  eV  for  electrons. 
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z  (cm) 

FIG.  7:  Benchmark  of  implicit  algorithm  vs.  explicit  algorithm.  Falsecolor  map  of  the  electron 
temperature  at  t  =  21  ns  for  (a)  explicit  algorithm  and  (b)  implicit  algorithm.  Note  that  the  radial 
scale  becomes  highly  nonlinear  for  cell  indices  greater  than  50,  as  illustrated  in  Fig.  [8} 

A.  Streamer  Propagation 

Effective  streamer  propagation  over  long  distances  relies  on  field  enhancement  which 
in  turn  depends  on  the  quality  of  the  streamer  as  a  conductor.  One  criterion  for  a  good 
conductor  is  that  the  electrostatic  shielding  distance  should  be  much  shorter  than  the  typical 
dimension  of  the  conductor.  In  our  modeling  we  indeed  found  that  streamers  propagated 
better  if  they  satisfied  this  criterion.  For  this  reason,  the  streamer  simulation  presented  here 
uses  a  radius  of  0.7  mm,  which  is  far  larger  than  the  radius  of  the  pre-ionization  associated 
with  a  femtosecond  laser  filament  [20] .  The  other  parameters  for  the  simulation  can  be 
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FIG.  8:  Nonlinear  mapping  between  radial  cell  index  and  position  for  use  in  interpreting  Fig.  [7] 


found  in  Table  |III|  The  simulation  ran  for  approximately  16  hours  on  64  processors  of  the 
Cray  XT3  “Sapphire”  at  the  Army  Engineer  Research  and  Development  Center. 

The  chemical  reactions  in  Table  III  are  taken  from  Refs.  [15]  and  [T9],  with  the  following 


caveats.  First,  the  excitation  rate  is  the  aggregate  rate  for  nitrogen  triplet  states.  The  actual 
de-excitation  mechanisms  are  not  modeled.  Instead,  a  fictitious  collisional  de-excitation 
mechanism  which  puts  the  excitation  energy  directly  into  the  gas  temperature  is  used. 
The  rate  associated  with  this  interaction  is  chosen  to  give  gas  heating  rates  similar  to 
those  observed  in  CHMAIR  simulations.  Second,  the  recombination  rates  technically  apply 
to  dissociative  recombination.  However,  since  dissociated  species  are  not  included  in  the 
model,  it  is  assumed  that  the  dissociated  species  immediately  re-associate,  and  that  the 
ionization  energy  goes  into  heating  the  gas.  Finally,  attachment  processes  have  not  been 
included  in  this  example.  This  can  be  partially  justified  by  the  fact  that  in  an  actual 
discharge,  photo-detachment  will  tend  to  offset  attachment.  In  fact,  visible  photons  are 
easily  energetic  enough  to  detach  an  electron  from  an  oxygen  anion  via  a  single  photon 
process  (the  attachment  energy  is  0.44  eV).  Since  a  real  streamer  emits  visible  photons, 
such  photo-detachment  is  indeed  expected  to  occur. 

The  on- axis  electric  field  is  plotted  as  a  function  of  z  for  various  times  in  Fig.  [9]  As 
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TABLE  III:  Parameters  for  Streamer  Propagation  Simulation 


Parameter 

Value 

Comment 

Electrode  Voltage 

500  kV 

constant  voltage 

Electrode  Diameter 

40  cm 

spherical  electrode 

Electron-Neutral  Cross  Section 

5  x  10~15  cm2 

ns,  Streamer  Density 

1014  cm-3 

see  Eq.  ( 

65 

) 

ro,  Streamer  Radius 

700  pm 

1/e  Definition 

Ls,  Streamer  Length 

1  cm 

Lt  —  Ls 

= 

.00  pm 

rif,  Preionization  Density 

1012  cm-3 

see  Eq.  ( 

65 

) 

iih,  Halo  Density 

0 

7ib,  Background  Density 

109  cm-3 

Ionization  Ratea 

2.2  x  10~8e_14-8/Te  cm3/s 

e~  +  02 

- 

2e"  + 

Ionization  Rate 

1.6  x  10“8Te^2e_17'2/Te  cm3/s 

e~  +  N2 

2e~  +  A2+ 

Excitation  Rate6 

5.4  x  io-7T“°'32e-9"52//Te  cm3/s 

e~  +  N2 

- 

e_  +  N2  [*] 

De-excitation  Rate 

3  x  10”9  cm3/s 

N2  [=t=]  +  n2 

— >  2N2 

Recombination  Rate 

0.0019T-°-5  cm3/s 

e-  +  0+ 

o2 

Recombination  Rate 

4.3  x  10”8Te_0'39  cm3/s 

e-  +  N+  - 

>n2 

Vibrational  Excitation 

see  Ref.  [15] 

Nz,  Axial  Cells 

3200 

Nr.  Radial  Cells 

150 

Nu,  Uniform  Cells 

75 

see  Eq.  ( 

63 

) 

At,  Time  Step 

3.55  ps 

“Temperatures  are  in  eV 

“7V2  [*]  represents  any  triplet  state 

expected,  strong  field  enhancement  produces  a  spike  in  the  electric  field  indicating  the  tip 
location.  At  71  ns,  this  held  enhancement  produces  a  held  well  in  excess  of  the  breakdown 
held.  Later  in  time,  the  held  enhancement  is  reduced  as  the  streamer  slows  down.  It 
continues  to  propagate  even  when  the  held  is  less  than  the  breakdown  held  due  to  the  fact 
that  attachment  was  neglected.  Fig.  [To] shows  the  electron,  vibrational,  and  gas  temperatures 


after  680  ns.  The  electron  temperature  tends  to  follow  the  electric  held,  as  is  expected  from 
theory  provided  the  vibrational  temperature  remains  low.  As  shown  in  the  figure,  the 
vibrational  temperature,  Tv,  is  indeed  much  lower  than  the  electron  temperature,  Te.  Thus, 
the  leader  phase  (defined  by  Tv  &  Te)  is  never  reached  in  this  simulation. 


FIG.  9:  Simulation  of  streamer  propagation:  on-axis  Ez  vs.  z  at  various  times. 


B.  Streamer  to  Leader  Transition 


To  model  the  streamer  to  leader  transition,  the  40  /irn  radius  typical  of  a  femtosecond  laser 
ionized  filament  is  used  as  the  pre-ionization  radius.  The  voltage  on  the  40  cm  diameter 
sphere  starts  at  500  kV  and  is  ramped  up  at  the  rate  of  570  V/ns.  Other  simulation 


parameters  are  given  in  Table  IV  Note  that  Coulomb  collisions  and  attachment  processes 
are  included. 

The  on-axis  electric  held  is  plotted  as  a  function  of  z  for  various  times  in  Fig.  [TT]  The 
evolution  of  the  spike  in  the  electric  held  is  markedly  different  from  the  case  of  the  streamer 
simulation  of  Fig.  [9]  The  held  enhancement  early  in  time  is  less  effective  because  the 
diameter  of  the  plasma  is  small,  and  therefore  there  is  less  charge  available  to  perturb 
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TABLE  IV:  Parameters  for  Streamer  to  Leader  Transition 


Parameter 
Electrode  Voltage 
Electrode  Diameter 
Electron-Neutral  Cross  Section 
Ion-Neutral  Cross  Section 
Electron-Ion  Cross  Section 
ns,  Streamer  Density 
r o,  Streamer  Radius 
Ls ,  Streamer  Length 
rif,  Preionization  Density 
rih,  Halo  Density 
rib,  Background  Density 
Ionization  Ratea 
Ionization  Rate 
Excitation  Rate6 
De-excitation  Rate 
Recombination  Rate 
Recombination  Rate 
Vibrational  Excitation 
Attachment  Rate 
Attachment  Rate 
Detachment  Rate 
Nz.  Axial  Cells 
Nr,  Radial  Cells 
Nu,  Uniform  Cells 
At,  Time  Step 


Value 


500  kV 


40  cm 


5  x  10”15  cm 
2  x  10~17 


cm 


see  Eq.  (18|) 
1016  cm”3 


2 

2 


40  /im 
1  cm 

1015  cm”3 
0 

109  cm”3 

2.2  x  iO”8e”14-8/Te  cm3/s 

1.6  X  10”8Te^2e”17'2/Te  cm3/s 

5.4  x  10”'Te”°'32e”9'52/Te  cm3/s 
3  x  10”9  cm3/s 
0.0019T”0-5  cm3/s 

4.3  x  10”8Te”°'39  cm3/s 
see  Ref.  [15i 

3.5  x  10”312r1e-°-052/Te  cm6/s 
10”31  cm6/s 

4.8  X  10 -WT1.5e-0A3/Tg  cm3/g 
800 
150 
50 


2.84  ps 


Comment 

570  V/ns  slew  rate 

spherical  electrode 


see  Eq.  (65 ) 


1/e  Definition 
Lf  —  Ls  =  100  fin r 


see  Eq.  (65 ) 


e”  +  02  — * > 

2e”  +  0^ 

e”  +  N2  — * ► 

2e”  +  1V2+ 

e”  +  N2  — >■ 

e”  +  2V2  [*] 

-Vjf*]  +  N2 

->  2  N2 

e”  +  0+  ^ 

■O2 

e”  +  N+  - 

>n2 

e  +  O2  +  O2  — *  O2  +  O2 

e”  +  O2  +  N2  — ^  N2  +  02 

N2  +  02  - 

>  O2  +  N2  +  e 

see  Eq.  (63) 


“Temperatures  are  in  eV 
bV2[*]  represents  any  triplet  state 
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FIG.  10:  Simulation  of  streamer  propagation:  on-axis  temperatures  vs.  2  after  680  nanoseconds. 


the  applied  field.  Later  in  time,  the  held  enhancement  increases  because  of  the  effects  of 
attachment.  Attachment  reduces  the  electron  density  ahead  of  the  streamer  which  leads  to 
a  more  extreme  gradient  in  charge  density  at  the  tip.  This  causes  the  electric  held  spike  to 
become  narrower  and  larger  in  amplitude. 


The  on-axis  electron,  vibrational,  and  gas  temperatures  after  220  ns  are  plotted  in  Fig.  12 


The  voltage  on  the  sphere  at  this  time  is  625  kV.  The  vibrational  and  gas  temperature  are 
approximately  1  eV  near  the  electrode.  The  vibrational  temperature  drops  to  0.5  eV  about 
8  cm  from  the  electrode.  The  streamer  tip  can  be  seen  as  the  spike  in  electron  temperature 
about  15  cm  from  the  electrode.  Since  the  vibrational  temperature  is  high  throughout  a 
significant  portion  of  the  filament  body,  the  transition  to  a  leader  can  be  considered  to  have 


occurred.  Figure  13  plots  the  tip  position  as  a  function  of  time.  For  the  hrst  100  ns  of 
propagation,  the  tip  slows  down  much  as  in  the  pure  streamer  simulation  discussed  above. 
After  100  ns,  the  tip  propagates  at  a  constant  speed,  and  even  appears  to  speed  up  near  the 
end  of  the  simulation.  This  behavior  will  be  investigated  further  as  this  work  continues. 

Although  this  result  is  encouraging,  it  does  not  illustrate  the  full  leader  development. 
The  missing  element  is  the  reduction  of  the  electric  held  in  the  leader  body.  The  electric 
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FIG.  11:  Simulation  of  streamer  to  leader  transition:  on-axis  Ez  vs.  z  at  various  times.  Note  that 
the  applied  field  increases  with  time. 


z  (cm) 


FIG.  12:  Simulation  of  streamer  to  leader  transition:  on-axis  temperatures  vs.  z  after  220  ns. 
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Time  (ns) 


FIG.  13:  Simulation  of  streamer  to  leader  transition:  position  of  streamer/leader  tip  vs.  time. 

field  in  the  body  is  a  constant  15  kV/cm,  which  is  too  high  to  allow  for  long  propagation 
distances.  This  may  be  due  to  the  small  diameter  of  the  leader  body  which  implies  a  large 
resistance.  It  may  be  that  including  hydrodynamic  effects,  or  running  the  simulation  longer, 
will  allow  the  leader  body  to  expand  radially  so  that  the  resistance  is  reduced. 

To  illustrate  the  hydrodynamics  capability  of  SPARC,  a  short  run  was  made  with  heavy 
particle  motion  enabled.  The  parameters  are  as  in  the  leader  simulation  above,  except  that 
electron-ion  collisions  are  neglected,  attachment  processes  are  neglected,  and  the  voltage 
on  the  sphere  is  held  fixed  at  +500  kV.  In  this  simulation,  heating  of  the  gas  results  in 
hydrodynamic  expansion  which  leads  to  a  noticeable  drop  in  the  on-axis  nitrogen  density 


within  70  ns.  This  is  illustrated  in  Fig.  14 


VII.  CONCLUSIONS 

A  simulation  code  has  been  developed  which  models  electrical  discharges  in  air.  The  code 
is  capable  of  modeling  streamer  propagation  over  long  distances  because  of  an  adaptive  grid, 
fully  scalable  parallelism,  and  an  implicit  model  for  electron  motion.  An  example  in  this 
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FIG.  14:  Simulation  with  heavy  particle  motion:  on-axis  nitrogen  density  vs.  z  after  70  nanosec¬ 
onds. 

report  shows  a  streamer  propagating  for  40  cm.  The  code  is  also  capable  of  modeling  the 
vibrational  excitation  of  nitrogen  and  its  interplay  with  ohmic  heating  of  the  electrons. 
This  allows  the  leader  phase  to  develop.  An  example  in  this  report  shows  the  vibrational 
temperature  reaching  about  1  eV  in  a  several  centimeter  long  region.  However,  a  dramatic 
reduction  in  the  on-axis  electric  held  was  not  observed.  Finally,  a  benchmark  comparing  the 
explicit  and  implicit  models  showed  that  the  implicit  model  recovers  the  important  features 
of  the  explicit  model  in  l/30th  of  the  computer  time. 
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Appendix:  TurboWAVE  Normalizations 

Although  this  report  is  written  in  cgs  units,  SPARC  is  not.  SPARC  consists  of  modules 
built  on  the  turboWAVE  framework.  This  framework  was  originally  designed  to  treat  colli- 
sionless  plasmas,  and  assumes  a  certain  normalization  scheme  appropriate  for  such  problems. 
Although  this  scheme  is  not  as  appropriate  for  a  code  like  SPARC,  it  is  adopted  for  the  sake 
of  keeping  everything  built  on  the  framework  consistent. 

The  unit  of  time  is  ujpl  where  ujp  =  2.8  x  1014  rad/s  is  the  plasma  frequency  in  fully 
ionized  air.  The  unit  of  length  is  c/up,  where  c  is  the  speed  of  light.  The  unit  of  mass  is  the 
electronic  mass,  m,  and  the  unit  of  charge  is  the  electronic  charge,  |e|.  The  unit  of  density  is 
np  =  mujp/ 47re2 .  Note  that  this  results  in  the  peculiarity  that  the  unit  of  particle  number  is 
Np  =  mc3/4:7iu}pe2.  To  determine  the  value  of  a  normalized  quantity  in  either  SI  or  gaussian 
units,  multiply  the  normalized  quantity  by  the  value  given  in  Table  [Vj  For  the  case  of  SI 
units,  we  define  the  permittivity  of  free  space,  eo  =  8.85  x  1CT12  F/m,  and  the  impedance 
of  free  space,  r/0  =  377  hi. 
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TABLE  V:  Normalization 


Quantity 

Symbol 

SI  Unit 

cgs  Unit 

Unit  Value 

Time 

(jp1 

(jp1 

3.55  fs 

Length 

c/ L Op 

c/ (Jp 

1.06  /mi 

Density 

np 

e^mUp/ e2 

mu2  /  47re2 

2.5  x  1019  cm“3 

Particle  Number 

eo  me3  /(jpe2 

me3 /Aiujpe2 

3  x  107 

Electric  Field 

Ep 

mcujp/e 

meup/e 

4.81  GV/crn 

Electric  Potential 

vie2 /e 

vie2 /e 

512  kV 

Magnetic  Field 

mojp/e 

majp/e 

1590  T 

Current  Density 

jp 

npec 

npec 

1.2  x  1011  A/cm2 

Charge  Density 

npe 

npe 

4  C/crn3 

Conductivity 

jp/  Up 

jp/  Ep 

2500  mho/m 

2.3  x  1013  s”1 

2-Body  Coefficient 

(Op/  Tip 

(Jp/ Tip 

1.1  x  10~5  cm3/s 

3-Body  Coefficient 

(jp/n2p 

(jp/n2p 

4.5  x  10~25  cm6/s 

Mobility 

c/Ep 

c/Ep 

6.2  cm2/V-s 

Energy 

me 2 

vie2 

8.2  x  10'14  J 

Energy  Density 

me2  np 

Vic2  Up 

2.1  MJ/crn3 

Power 

mc2iop 

VlC2(Jp 

23  W 

Power  Density 

vic2(jpnp 

VlC2(JpVp 

5.77  x  102°  W/crn3 

Temperature 

vie2 /ks 

vie2  /ks 

512  keV 

5.93  x  109  K 

Fluence 

El/m^p 

cEp/4irtjp 

220  J/crn2 

Intensity 

Ip 

Elho 

cE2  /  47t 

6.1  x  1016  W/crn2 

Cross  Section 

L Op  j TlpC 

C Op  j TtpC 

3.76  x  10~16  cm2 

